Geographic variation in evolutionary rescue under climate change in a crop pest–predator system

Abstract Species distribution models (SDMs) are often built upon the “niche conservatism” assumption, such that they ignore the possibility of “evolutionary rescue” and may underestimate species' future range limits under climate change. We select aphids and ladybirds as model species and develop an eco‐evolutionary model to explore evolutionary rescue in a predator–prey system under climate change. We model the adaptive change of species' thermal performances, accounting for biotic interactions. Our study suggests that, without considering evolutionary adaptation, the warming climate will result in a reduction in aphid populations and the extinction of ladybirds in large parts of the United States. However, when incorporating evolutionary adaptation into the model, aphids can adapt to climate change, whereas ladybirds demonstrate geographic variation in their evolutionary rescue potential. Specifically, ladybirds in southern regions are more likely to be rescued than those in the north. In certain northern regions, ladybirds do not avoid extinction due to severe warming trends and seasonality of the climate. While higher warming trends do prompt stronger evolutionary changes in phenotype, they also lead to reduced aphid population abundance such that ecology constrains ladybird population growth. Higher seasonality induces an ecological effect by limiting the length of reproductive season, thereby reducing the capacity for evolutionary rescue. Together, these findings reveal the complex interplay between ecological and evolutionary dynamics in the context of evolutionary adaptation to climate change.

Eq. 6 The birth rate for an aphid individual at time t .Eq. 6 The mean birth rate for the aphid population at time t .Eq. 6 The intrinsic mortality rates for an aphid individual or ladybird individual at time t , respectively.
Eq. 6 The mean intrinsic mortality rates for aphids or ladybirds at time t , respectively.Eq. 6 The effect of temperature on the functional response, range from 0-1.Eq. 8 The instantaneous rate of growth of the descendant lineage for an aphid with phenotype z A,i at time t .
Eq. 9 The instantaneous rate of growth of the descendant lineage for a ladybird with phenotype z L,j at time t .
Eq. 11 The probability densities of aphid or ladybird individuals with genotype z g ,A or z g ,L .
Eqs. 21-22 The probability densities of aphid or ladybird individuals with environmental effect z e,A or z e,L .

Eqs. 21-22 ρ τ
The probability density for the mutational effect of aphids.Eq. 30 Eq. S1 CT min , CT max Minimum or maximum temperature thresholds beyond which growth rate is nil, mortality rate is maximal.CT min = z A,i /z L,j − 15, CT max = z A,i /z L,j + 10.

Eq. S1
CT opt1 ,CT opt2 Temperature range in which mortality is minimal.CT opt1 = z A,i /z L,j − 5, CT opt2 = z A,i /z L,j + 5 Eq. S3 q 1 , q 2 Shape parameters used to adjust the skewness of the thermal performance curve.
q 1 = 1.5, q 2 = 1 Eq. S1 Eq.S2 perature range in which mortality is minimal.q 1 and q 2 are shape parameters which adjust the skewness of the 18 thermal performance curve.k 1 , b 1 , k 2 and b 2 are the parameters which make υ (CT opt1 ) = υ (CT opt2 ) = υ max and 19 υ (CT min ) = υ (CT max ) = υ mi n .We hold q 1 = 1.5 and q 2 = 1 constant, and set CT min , CT max , CT opt1 and CT opt2 as 20 functions of T opt , to ensure the thermal performance curves have the same shape over time.For g (z A,i , T , t ) and in 21 the context of fecundity, m is set to be m A = 0.6 to roughly match the maximum intrinsic rate of growth of aphid

23
The impact of temperature on the predation rate of ladybirds is already accounted for in Eq. 8, which uses 24 g (z L,j , T , t ) to adjust the functional response of ladybirds, setting m L = 1.In addition, the temperature effect on 25 the birth rate of ladybirds is represented through the transformation rate Q p (i.e., the mean number of aphids a lady-26 bird needs to consume to reproduce a single egg), which determines the numerical response of ladybirds.Thus, the 10 temperature-dependent Q p for each individual ladybird with phenotype z L,j is  For the aphid population, the un-normalized probability density of each aphid genotype is To model the population dynamics of aphids and ladybirds in various geographic locations, we used daily temperature 40 data in 2000 to estimate the locally adapted optimal temperature for each species ( z * A for aphids and z * L for ladybirds).

41
These estimates are used as the starting values for our long-term simulations.We followed the procedures below to 42 estimate z * A and z * L .

43
To begin, we obtained the maximum daily temperature data in 2000 (T max,2000 ) and used its floor value (i.e., the 44 greatest integer that is less than or equal to T max,2000 ) as the initial values for z * A and z * L .
zA and zL .If the trend value is positive, it means that the value of z * A is too low and needs to be increased.Conversely, a negative trend value indicates that z * A is too high and needs to be decreased.To determine the optimal value of z * A , T opt for aphid and ladybird  !"#,%  !"#,& F I G U R E S 1 Changing trend of annual averaged T opt for aphid (T opt,A , blue curve) and ladybird T opt,L , orange curve) over 150 years.d) Evolutionary potential for ladybirds ( ,, ) g) Evolutionary potential for aphids ( ,, ) Daily aphid population abundance and effect of evolutionary potential on aphid and ladybird population abundance (AAP and ALP) and their evolving trait, T opt , in the 'rescued' (SO: Logansport, Louisiana) and 'extinct' (NO: Beauford, Minnesota) locations.(a) daily population abundance for aphids from 2050-2100 in the two locations.(b1)-(d2) represent the changing trend for AAP, ALP and T opt when altering ladybirds' segregration variance (Black curves:V g ,L,0 =0.3; Red curves:V g ,L,0 =3).(e1-g2) represent the changing trend for AAP, ALP and T opt when altering aphids' segegration variance (Blue curves: V g ,A,0 = 0.15, Black curves:V g ,A,0 =0.3; Red curves:V g ,A,0 =3) for ladybirds.In (d1), (d2), (g1) and (g2), solid curves represent T opt,A , dashed curves represent T opt,L .U R E S 3 Ladybird population abundance (ALP) in 2150 under different climates with various seasonality (s) and warming trend (k ).(a)-(c) represent three different locations (SO: Logansport, Louisiana; CO: Taberville, Missouri; NO: Beauford, Minnesota).The color depth in the heatmap indicates the magnitude of ALP, which also represents the degree of evolutionary rescue.A lighter color indicates weaker evolutionary rescue.The black box in each subplot indicates the actual combination of seasonality and trend in each location.The orange numbers in each subplot denote the rate of reduction in ladybird population abundance as we increase seasonality or trend.

T 1 s 2 k
Daily temperatures and ladybird density in 2000 and 2100 in Louisiana (SO).'Baseline' (blue curves) means we use the actual seasonality value in Louisiana.'Higher seasonality' (orange curves) represents that we use the seasonality value in Minnesota (NO), which has a much higher seasonality.(a)-(b) daily temperature in 2000 and 2100.(c)-(d) daily ladybird density in 2000 and 2100.We only plot the reproductive season in (c) and (d).TA B L E S 1 Notations and parameters in the eco-evolutionary model.Initial annual mean temperature in 2000 in each location.Eq.Seasonality, defined as half of the difference between yearly minimum temperature T mi n,y e ar and yearly maximum temperature T max,y e ar Eq.Warming trend, defined as the change of annual mean temperature every century in each location Eq. 3 DL Daylength, the required number of daylight hours, is determined by ξ, which is the exposed radius between the sun's zenith and the sun's solar circle and estimated by a function of latitude (φ), the number of days from January 1st (t ) and the Earth's rotational axis R = 23.439• Eqs.36 ε A computational control "switch" to enable the sexual phase, aphids remain in their asexual phase when ε = 0 and switch to their sexual phase when ε = 1.Eq. 36 A, L Population density of aphids or ladybirds -AAP, ALP Accumulation of the daily population abundance for aphids or ladybirds within each year z A,i , z L,j Phenotypic value of optimumal temperature T opt for an aphid individual or ladybird individual.Eqs.4-5 z g ,A,i , z g ,L,j Genotypic value of an aphid individual or ladybird individual Eqs.4-5 z e,A,i , z e,L,j Environmental effects on the phenotype of an aphid individual or ladybird indi-value Note zg,A , zg,L The mean breeding value for the population of aphids or ladybirds.-zA , zL The mean phenotypes for the population of aphids or ladybirds z * A , z * L Locally adapted optimal temperature for aphids or ladybird -{z A,i }, {z L,j } Continuous set of phenotypes in the population for aphids or ladybirds, respectively.

η
genetic variance for aphids or ladybirds.Segregation variance V g ,A,0 = V g ,L,0 = 0.3 -Continued on next pageTable S1 -continued from previous page Notation Definition and Parameter value Note V e,A , V e,L Environmental variances for aphids and ladybirds.V e,A = V e,L = 0.7 f Inbreeding coefficient for the aphid population.f = 0Mutational heritability for aphids (s = A) and ladybirds (s = L).We selected 0.002 (estimated for the thermal performance of Drosophila) as the approximate estimate for h 2 m,A and h 2 m,L (Houle et al., 1996; Lynch and Walsh, 1998; Latimer et al., 2014).-Generation time coefficient for aphids (s = A) and ladybirds (s = L).The rough estimates for η A and η L are 0.05 generations/half-day and 0.01 generations/halfday, respectively.-V τ,s Mutational variance for aphids s = A and ladybirds (s = L).V τ,s = h 2 m,s × V e,s × η s .Eqs. 25, 34 K Carrying capacity for aphid population.Eq. 6 m A Optimal growth rate at optimum temperature, its value depends on temperaturedependent generation time.m A = 0.6
S3) where Q * p = 500 and Q ′ p = 2000 represent the minimum and maximum values for Q p which ensures the range of 29 the intrinsic rate of growth for ladybirds (r mL : 0 − 0.2) matches with empirical data under ideal condition (i.e., prey 30 saturation and no competition) (Islam et al., 2022).In our case, r mL is typically in the range of 0 − 0.02 given the 31 carrying capacity of aphids in our model (Fig. S5).

32
Aphid population density (A) Intrinsic rate of growth for ladybirds (!" ) Theoretical maximum  !" ≈ 0.2 Carrying capacity ( = 5×10 ! ) Within the range of  F I G U R E S 5 Intrinsic rate of growth for ladybirds changes with aphid population density.

S1. Parameters for temperature-dependent vital rates.
In our model, all of the thermal traits (e.g., CT min and CT max ) for aphids and ladybirds are assumed to evolve with 5 climate change at the same rate with T opt , so they can be expressed as functions of z A,i or z L,j and time t .z A,i T is the air temperature.CT min and CT max are the minimum and maximum temperature thresholds 16 beyond which birth rate (or predation rate) is nil, or mortality rate is maximum.CT opt1 and CT opt2 define the tem- 4 8 temperature-dependent functions within the stage-structured population dynamic model previously developed by Ge 9 et al. (2023) to estimate the vital rate at time t for individual aphid or ladybird.The functions for an aphid individualg (z A,i , T , t) is used to denote the temperature-dependent birth rate for an aphid individual at time t , g (z L,j , T , t ) de-13 notes temperature effect on the predation rate and birth rate of a ladybird individual at time t .In Eq.S2, υ (z A,i , T , t ) 14 and υ (z L,j , T , t ) are utilized to represent the intrinsic mortality rate for an aphid individual and ladybird individual, 15 respectively.